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ABSTRACT 

An  analytic  model  is  developed  of  baroclinic  currents 
measured  by  several  investigators  under  the  Arctic  Ocean  ice 
cover.   A  time-dependent,  Boussinesq,  quasi-geostrophic 
vorticity  equation  is  developed  to  describe  the  flow. 
Separable  normal  mode  solutions  are  found  which  show  sta- 
tionary current  patterns  qualitatively  consistent  with 
observed  data.   An  initial  value  problem  is  then  formulated 
using  density  perturbations  at  a  newly  forming  lead  as  a 
forcing  function. 
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I.   INTRODUCTION 

A.   OBJECTIVES  OF  RESEARCH  AND  REVIEW  OF  PREVIOUS  RFSFARCH 
The  objective  of  this  study  was  to  develop  an  analytic 
model  of  specific  small  scale  baroclinic  current  patterns 
which  have  been  observed  on  a  recurring  basis  beneath  the 
Arctic  ice  cap.   The  patterns  are  termed  small  in  that  they 
are  on  a  scale  less  than  that  of  the  general  Arctic  circula- 
tion.  In  order  to  develop  the  model  a  theoretical  investi- 
gation was  made  of  the  baroclinic  response  of  the  Arctic 
Ocean  to  perturbations  in  its  internal  density  structure  at 
a  newly  forming  lead. 

It  is  worthwhile  at  this  point  to  briefly  discuss  pre- 
vious research  in  the  area  of  Arctic  dynamic  modeling,  and 
to  describe  the  baroclinic  currents  measured  by  several 
investigators,  in  an  effort  to  indicate  the  relationship  of 
this  study  to  previous  work. 

Investigations  of  Arctic  Ocean  dynamics  can  be  classified 
by  their  scale  as  well  as  by  the  mode  of  motion  described 
(either  internal,  surface,  or  a  combined  mode).   The  largest 
scale  studied  concerns  processes  which  affect  the  entire 
ocean,  such  as  long  term  mean  ice  drift  and  water  circulation 
In  order  to  practically  study  motion  on  this  scale,  spacing 
between  grid  points  cannot  feasibly  be  less  than  tens  of 
kilometers.   This  necessitates  parameterizing  some  of  the 
smaller  scale  processes  occurring  between  grid  points,  rather 
than  solving  for  them  explicitly.   Factors  which  must  be 
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considered  on  this  scale  are  the  shape  and  bathymetry  of  the 
ocean  basin,  Coriolis  effects,  wind  and  ice  stress,  etc. 

Several  theoretical  investigations  of  large  scale  Arctic 
Ocean  dynamics  which  closely  match  observed  patterns  of  water 
circulation  and  ice  drift  have  been  carried  out  in  recent 
years.   Among  these,  Campbell  [1965]  proposed  a  large  scale 
analytic  ice  drift  model  and  Gait  [ 197 3a]  developed  a  working 
barotropic  numerical  model  of  overall  surface  water  circula- 
tion.  Both  these  models  were  concerned  with  surface  motion 
or  an  equivalent  barotropic  circulation.   A  complete  model 
of  large  scale  internal  /arctic  circulation  has  not  yet  been 
published,  but  work  along  these  lines  is  being  carried  out 
at  NOAA's  Geophysical  Fluid  Dynamics  Laboratory,  Princeton, 
New  Jersey  [Semtner,  1973]. 

The  smallest  scale  of  analysis  concerns  the  short  term 
local  effects  of  phenomena  such  as  refreezing  leads  and  ice 
cracking.   These  effects  result  in  dynamic  processes  such  as 
forced  convection  due  to  the  salt  rejection  and  heat  transfer 
during  freezing.   Grid  size  of  concern  here  would  be  from 
less  than  a  cm,  a  microscale,  to  a  few  meters.   Coriolis 
effects  would  not  be  important  at  the  smallest  end  of  this 
scale  but  effects  such  as  diffusion  would  have  to  be  con- 
sidered.  Several  investigators  have  done  work  on  the  problem 
of  forced  convection  associated  with  freezing  of  sea  ice. 

Foster  [1972]  treated  the  problem  analytically  after 
having  published  several  earlier  papers  on  the  problem  of 
forced  convection  in  the  ocean.   The  initial  state  in  his 


study  was  a  homogeneous  under-ice  layer.   He  then  assumed 
density  to  be  a  function  of  salinity  only  and  solved  the 
equations  of  motion  and  salt  diffusion  directly  associated 
with  a  refreezing  lead.   Using  a  two-dimensional  model  Foster 
carried  his  solution  to  the  point  of  finding  initial  vertical 
velocities  and  length  and  depth  scales  for  convection  cells 
during  the  initial  stages  of  convection.   These  were  about 
1.4  cm/sec  for  vertical  velocities,  80  cm  for  horizontal 
convection  cell  spacing,  and  several  cms  for  depth  of  con- 
vective  penetration.   He  did  not  analytically  consider 
further  stages  of  convection  nor  did  he  continue  to  a  larger 
scale  in  which  Coriolis  forces  would  become  important.   He 
did  postulate  a  higher  order  of  convection  cells  with  depth 
penetration  of  up  to  6  meters  and  vertical  velocities  of 
about  1  cm/sec. 

Piascek  [1970]  considered  the  general  problem  of  surface 
cooling  effects  on  the  ocean.   He  developed  a  numerical 
model  for  convection  cells  on  both  a  small  scale  (order  of 
cms)  which  neglected  Coriolis  effects  and  a  larger  scale 
(order  of  kms)  which  considered  Coriolis.   His  model  con- 
figuration was  a  rectangular  box  directly  below  the  surface 
with  an  initially  homogeneous  structure.   He  considered 
uniform  surface  cooling  as  a  boundary  condition  and  then 
developed  a  vorticity  equation  in  terms  of  a  vector  poten- 
tial.  This  was  numerically  integrated  over  a  2-dimensional 
grid  to  give  steady  state  temperature  and  vorticity  patterns 
within  his  "box."   With  some  modification  his  model  could  be 
adapted  to  describe  Arctic  under-ice  convection. 
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Coachman  [1966]  considered  a  mass  and  energy  balance  on 
the  layer  of  water  immediately  below  the  ice  cover.   He 
assumed  a  downward  flow  direction  under  the  refreezing  leads 
due  to  sinking  of  highly  saline  water  and  a  compensating 
upward  flow  spread  out  over  the  area  between  leads,  as  shown 
in  Figure  1.   His  study  used  mean  values  of  freezing  rate, 
temperatures,  lead  size  and  distance • between  leads  to  deter- 
mine mean  vertical  velocity  values.   His  scale  of  analysis 
was  necessarily  larger  than  that  of  Foster  in  order  to  arrive 
at  mean  velocities.   His  downward  vertical  velocity  values 
are  on  the  order  of  .01  to  .08  cm/sec.   The  depth  of  his 
surface  circulation  layer  was  100  meters.   If  the  Arctic 
Ocean  is  viewed  from  a  large  scale  standpoint,  the  net 
result  of  the  investigations  by  Coachman,  Foster,  and  Piascek 
would  be  to  show  that  there  is  a  small  downward  vertical 
velocity  associated  with  each  refreezing  lead  due  to  forced 
convection  and  that  there  is  a  local  increase  in  the  density 
of  the  water  directly  under  a  refreezing  lead. 

A  related  small  scale  problem  which  has  been  investigated 
is  that  of  the  ice  formation  rate  and  rates  of  salt  rejection 
and  heat  transfer  at  a  refreezing  Arctic  lead.  Schaus  [1971] 
developed  a  numerical  model  of  a  refreezing  lead  which  treats 
these  problems  and  which  could  be  used  in  coupling  the  air- 
sea  interaction  problem  with  ice  dynamics  and  with  the  under- 
ice  convection  problem  discussed  here. 

A  vertical  velocity  similar  to  that  caused  by  the  convec- 
tive  process  may  also  arise  from  the  initial  cracking  and 
ice  motion  during  lead  formation.   This  specific  problem 
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Figure  1.   Schematic  of  Convection  Model  Proposed 
by  Coachman  [19  6  6]  . 
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has  not  been  previously  studied.   In  this  case  the  theory  is 
as  follows,  assuming  ice  shears  upon  cracking  over  a  period 
of  a  few  days  as  shown  in  Figure  2.   There  is  a  stress 
imparted  to  the  underlying  water  at  the  ice-water  boundary, 
which  in  turn  causes  a  convergence  of  the  Fkrnan  layer.   This 
in  turn  results  in  an  impulsive  vertical  velocity  under  the 
newly  forming  lead.   The  net  result  is  that,  for  an  entirely 
different  reason,  a  vertical  velocity  similar  to  that  caused 
by  convection  is  associated  with  shear  across  newly  forming 
leads.   Both  these  lead-associated  phenomena  would  be  classed 
as  small  scale  effects. 

By  interacting  with  the  existing  density  and  flow 
patterns   these  velocity  and  density  perturbations  may  have 
some,  as  yet  undetermined,  effect  on  the  large  scale  internal 
flow  dynamics  of  the  Arctic.   A  proposed  theory  for  this  is 
discussed  and  developed  in  part  B  of  this  section.   Any 
resultant  flow  pattern  due  to  such  an  interaction  might 
then  be  expected  on  a  scale  intermediate  in  size  between 
the  smallest  and  largest  scales  discussed;  on  the  order  of 
meters  in  depth  and  meters  to  kilometers  in  length.   A  theo- 
retical investigation  of  baroclinic  internal  flow  in  the 
Arctic  on  this  scale  has  not  previously  appeared  in  the 
literature.   It  is  important  to  note  here  that  this  is  the 
scale  of  a  number  of  internal  current  patterns  or  events 
measured  in  the  Arctic. 

The  observed  currents  appear  to  be  baroclinic  in  nature 
although  a  satisfactory  theory  for  their  origin  has  not  yet 
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been  developed.   An  example  of  the  vertical  distribution  of 
horizontal  velocities  measured  by  Gait  [1967]  under  Ice 
Island  T-3  is  shown  in  Figure  3.   The  corresponding  density 
distribution  is  also  shown.   Similar  results  found  by 
Coachman,  and  Newton  [1973]  are  shovTn  in  Figure  4.   It  can  be 
seen  that  a  peak  in  the  horizontal  velocity  occurs  in  the 
pycnocline,  with  a  maximum  at  about  150  meters,  indicating 
a  baroclinic  origin.   Figure  5  shows  the  horizontal  distri- 
bution of  velocities  found  by  Coachman  and  Newton  at  150 
meters.   It  is  a  plot  of  velocity  vectors  measured  from 
one  position  on  the  ice  while  drifting  over  a  given  area. 
As  can  be  seen,  there  is  a  periodic  change  in  the  velocity 
direction  with  an  apparent  wavelength  of  about  30  km.   The 
reason  for  classing  this  periodicity  as  spatial  rather  than 
temporal  is  that  the  observed  currents  were  found  to  be 
essentially  in  geostrophic  balance.   For  this  type  feature, 
phase  propagation  speeds  will  be  small  compared  to  advective 
speeds.   Figure  6  shows  a  vertical  cross-section  of  veloci- 
ties as  reported  by  Coachman  and  Newton  [1973], 

A  need  thus  exists  for  an  investigation  of  the  Arctic 
Ocean  internal  flow  dynamics  on  the  intermediate  scale 
described.   The  immediate  purposes  of  this  study  are  to  offer 
an  explanation  of  the  observed  baroclinic  currents.   A 
further  use  of  such  a  study  is  to  couple  the  small  scale 
phenomena  to  the  overall  circulation,  which  in  turn  may  aid 
in  thfi  development  of  a  three  dimensional  Arctic  circulation 
model . 
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Figure  3.   Vertical  Distribution  of    Horizontal 
Velocities  Measured  bv  Gait  [1967]. 
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Figure  6.   Vertical  Section  of  Velocities  Measured 
by  Coachman  and  Newton  [19  73]  . 
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B.   THEORY  DEVELOPED  IN  THIS  THESIS 

A  theory  was  postulated  to  describe  the  dynamic  effects 
on  the  underlying  ocean  of  a  nev/ly  forming  Arctic  lead.   The 
theory  developed  was  begun  with  the  premise  that  a  small 
vertical  velocity  perturbation  or  density  increase  is  created 
during  the  formation  of  each  lead,  caused  by  either  refreez- 
ing  or  cracking,  as  described  earlier.   The  scale  of  analysis 
used  v/as  tens  of  kms  in  the  horizontal  and  the  full  depth 
of  water  (0-3000  m)  in  the  vertical. 

The  water  directly  below  the  Arctic  ice  cover  can  typi- 
cally be  expected  to  consist  of  a  shallow  homogeneous  layer 
(20-50  m) ,  under  which  a  stably  stratified  density  distribu- 
tion exists.   There  is  also  some  horizontal  slope  to  the 
isopycnals  in  order  to  support  a  mean  flow.   The  horizontal 
gradient  will  be  several  orders  of  magnitude  less  than  that 
of  the  vertical  density  gradient  in  the  pycnocline.   Pertur- 
bations such  as  the  vertical  velocity  and  density  changes 
occurring  under  a  lead,  which  is  a  very  small  area  on  the 
assumed  length  scale,  should  distort  the  isopycnals  over  a 
relatively  short  period  of  time  and  short  distance.   This 
would  act  as  an  impulsive  event  on  an  overall  time  scale 
measured  in  days.   As  the  isopycnals  attempt  to  return  to 
their  stable  condition,  a  transient  waveform  of  some  sort 
would  then  appear  and  propagate  out  from  the  lead. 

The  first  problem  would  then  be  to  find  the  form  of  the 
waves  which  can  exist  in  a  typical  Arctic  density  structure 
and  to  then  determine  if  these  are  consistent  with  observed 


20 


current  patterns.   Once  this  was  done,  a  second  problem  would 
be  an  investigation  to  determine  whether  the  perturbations 
at  a  lead  could  actually  force  these  oscillations.   An  inves- 
tigation of  these  two  problems  forms  the  basis  of  this  thesis. 

An  indication  of  the  possible  validity  of  the  work  at  any 
point  is  the  degree  of  correlation  of  theoretically  developed 
waves  or  currents  with  real  data.  A  comparison  was  therefore 
made  with  the  baroclinic  currents  measured  by  Gait  [1967]  and 
Coachman  and  Newton  [1973]  once  a  theoretical  current  pattern 
was  developed. 

The  first  step  in  substantiating  the  theory  was  to  inves- 
tigate the  type  of  waves  which  could  exist  in  the  given 
density  structure  and  to  compare  these  with  known  data.   This 
corresponded  to  treating  the  Arctic  Ocean  as  a  semi-enclosed 
system,  subject  to  certain  governing  equations  of  motion  and 
boundary  conditions ,  and  solving  for  its  preferred  modes  of 
oscillation.   The  actual  solution  to  this  problem  is  the 
subject  of  Section  II  of  this  thesis.   The  general  assumptions 
made  in  formulating  the  theory  and  the  rationale  behind  them 
will  be  discussed  here. 

The  equations  used  follow  the  development  of  Gait  [1973b] 
which  is  based  on  that  of  Kuo  [1973],   The  basic  difference 
from  Kuo's  development  is  in  the  formulation  of  the  equation 
of  state.   The  system  was  assumed  to  be  inviscid,  hydrostatic, 
and  quasi-geostrophic .   The  equation  of  state  assumes  density 
proportional  to  salinity  only.   This  is  a  reasonable  approx- 
imation in  the  region  being  investigated  as  can  be  seen  in 
Figure  7,  which  is  a  plot  of  salinity,  temperature,  and 
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Figure  7.   Observations  of  Salinity,  Temperature  and 
Density  Taken  Under  Ice  Island  T-3. 
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density  profiles  taken  under  T-3  along  with  current  measure- 
ments [Gait,  1967]  . 

The  Boussinesq  approximation  was  also  made,  which  states 
that  density  variations  are  significant  only  when  the  density 
is  multiplied  by  gravity.   An  additional  assumption  necessary 
to  this  development  is  that  the  Rossby  number  is  small.   This 
allowed  the  expansion  of  the  dependent  variables  as  power 
series  in  terms  of  the  Rossby  number.   Equations  consistent 
in  powers  of  Rossby  number  could  then  be  derived. 

From  the  above  assumptions  and  the  basic  equations  of 
motion  a  quasi-geostrophic  vorticity  equation  was  developed 
in  terms  of  the  perturbation  pressure  (pressure  could  be 
written  as  the  sum  of  a  mean  and  a  perturbation  quantity) . 
This  equation  was  solved  by  assuming  a  separable  solution 
consisting  of  the  product  of  a  horizontal  plane  wave  and  an 
unknown  vertical  structure  function.   This  vertical  structure 
function  was  then  solved  for  numerically  subject  to  appro- 
priate boundary  conditions.   Results  of  the  numerical  solu- 
tion are  compared  with  the  measured  currents  in  Section  III 
of  this  thesis. 

Investigation  of  the  second  problem  described  (whether 
the  theoretical  oscillations  could  actually  be  driven  by 
perturbations  found  at  a  lead)  was  accomplished  by  using 
pressure  conditions  under  a  newly  forming  lead  as  an  initial 
condition  for  the  quasi-geostrophic  vorticity  equation 
developed  in  Section  II.   A  description  of  this  initial  value 
problem  forms  Section  IV  of  this  thesis. 
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II.   DEVELOPMFNT  OF  THE  WAVE  SOLUTION 

A.   ANALYTIC  DEVELOPMENT  OF  THE  VORTICITY  EQUATION 

The  initial  formulation  of  the  problem  was  carried  out 
using  a  rectangular  coordinate  system,  with  positive  values 
of  x  to  the  East,  y  to  the  North,  and  z  upward.   The  govern- 
ing equations  of  motion  for  the  flow  can  be  written  in  the 
following  Boussinesq,  inviscid  form: 

u.  +  uu  +  vu   +  wu   -  f v  =  -  —  p  (1) 

t     x     y      z  p   x  x 

1  s 

v,  +  uv   +  w   +  wv   +  f  u  =  -  —  p  (2) 

t      x      y      z  ps   y 

w,  +  uw   +  vw   +  ww   =  -  —  p   -  —  g  (3) 

t     x     y      z      p   *z    p   ^ 

2  s        s 

where  u,v,  and  w  are  the  x,y,  and  z  velocity  components, 
respectively,  f  is  the  vertical  component  of  the  Coriolis 
parameter,  g  is  the  acceleration  of  gravity,  p  is  the  pres- 
sure, p  is  the  density  and  p   is  the  horizontal  space  and 
time  average  of  the  density,  which  is  a  function  of  z  only. 
In  addition  to  the  equations  of  motion,  the  continuity  equa^ 
tion  can  be  written  as: 

p.+up   +vp   +  wp   +p(u+v+w)=0       (4) 
t      x      y      z     s   x    y    z 

and  the  equation  for  conservation  of  salt,  assuming  a  con- 
servative lav;  with  no  mixing,  is: 

s.  +  us   +  vs   +  ws  =  0  (5) 

t      x      y      z 
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Equations  (1)  through  (5)  can  then  be  scaled  and  non- 
dimensionalized  using  essentially  the  same  system  as  was 
used  by  Kuo  [1973].   The  non-dimensional  quantities  (primed) 
are  defined  as  follows: 

(x,y)  =  L(x',y');         (u,v)  =  U(u',v');        (6) 

z  =  Dz  ;  w  =   — - —  w  ; 


\  L  I 


t 


=  <L-)t,;  f  =  fof'; 


P  =  Ps  +  f0ULps(}),  ;       p  =  ps  (1  +  ep'  )  ; 

s  =  s  ( 1  +  e  s  '  )  ; 

s 

where  L  is  a  characteristic  horizontal  scale,  D  is  a  charac- 
teristic depth  scale,  U  is  a  characteristic  horizontal 
velocity  and  p   and  s   are  the  static  pressure  and  density 
fields  respectively. 

In  addition  the  following  non-dimensional  coefficients 
are  defined: 

R  =  ■£—=-    ,   Rossby  number;  (7) 

0  J 

2  2 

f  ZL 

u  =  —  ,  a  ratio  of  the  scale  length,  L ,  to 

the  Rossby  deformation  length; 

£  =  uR  ,     a  measure  of  the  perturbations  in 
the  thermodynamic  quantities; 

1  8ps 

6   =  —  -r ,  a  measure  of  the  stability  of  the 

s      static  density  field; 

,   8s 
6   =  —  7; ,  a  measure  of  the  stability  of  the 

S     S    a  Z 

s      static  salinity  field. 
Substitution  from  (6)  into  (1)  and  (2)  then  gives: 
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R  (u  '  .  +  U  *  U  '  .  +  v  '  u  '  ,  +  Rw  '  u  '  ,  ) 


-  f  V'   =  -<|) 


X' 


(8) 


R  (v  '  ,  +  u  '  v  '  ,  +  v '  v  '  ,  +  Rw  '  v  *  )  +  f  '  u  ' 
t       x       y        z 


and  substitution  from  (7)  into  (4)  gives 


R2  f  ?.)   (w '  +  u  '  w  '  ,  +  v '  w '  ,  +  Rw  '  w  '  ,  ) 
\L/    t '      x '      y '       z 


y 


(9) 


z  '    pY 


providing  that  the  hydrostatic  relationship  for  the  steady 
state  density  and  pressure  fields  is  used,  as  given  by: 


(pj  =  -pQg 

S  z       s 


Then,  substitution  of  (6)  into  (4)  gives 


e(p£,  +  u'p^,  +  v'p^,  +  Rw'p'z,)  +  —  (w'ps)z,  +u^,  +  v^,  =  0  (11) 


Finally,  substitution  from  (6)    into  (5)  gives: 


(6 


R  (s  '  +  u  '  s  '  .  +  v  '  s  '  ,  +  Rw '  s  '   + 
t      x '      y        z 


M 


w '  )  -  0 


(12) 


Equations  (8)  through  (12)  can  then  be  simplified  by  consi- 
dering limitations  which  occur  in  the  real  ocean: 


D 


:<  1       used  in  eq.  (10) 


(13) 


e  =  yR  «  R  used  in  eq.  (11) 


(14) 


and  by  assuming  a  simplified  equation  of  state  relating 
density  to  salinity: 

p  =  s 

or,  in  a  non-dimensional  form  for  the  perturbation  quantities 


p'  =  s' ; 


6   =5   -  0  (R)  . 

P     s 


(15) 


2f 


Using  eq .  (15),  the  non-dimensional  equations  become 
(omitting  primes  from  here  on) : 

R(u,_  +  uu   +  vu   +  Pwu  )  -  fv  =  -<J>  (16) 

t      x      y       z  ^x 

R(v^  +  uv   +  vv   +  Rwv  )  +  fu  =  -4  (17) 

t      x      y       z  y 

<j>   +  6  6  +  s  =  0  (18) 

Yz     pY 

u   +  v   +  —  (wp  )   =0  (19) 

x    y   ps    s'z 

R(s.  +  us   +  vs   +  Rv;s   +   -2-  w)  =  0  (20) 

t      x      y       z    1  \i  I 

At  this  point  the  assumption  of  a  small  Rossby  number,  R, 
will  be  used  and  dependent  variables  expanded  in  a  power 
series  form;  i.e.,  u,  v,  6,  and  s  will  take  the  form: 

a  =  E   Rna  (21) 

n=0     n 

where  the  a  represents  the  general  dependent  variable  and  w 
will  be  written  as: 

2 
w  =  w,  +  Rw„  +  R'w-,  +  ...  etc. 

This  is  consistent  with  the  scaling  given  by  equation  (6). 
In  addition,  the  Coriolis  parameter  will  be  scaled  on  a 
beta-plane : 

fQ  =  1  +  By  (22) 

where 

By  <   0(R)  (23) 

Using  these  approximations  plus  eg.  (15)  and  the  series 
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expansion  gives  the  following  equations  to  zero  order  in 
Rossby  number: 

-v   =  -(J)  (24) 

o     o 

X 

u  =  -4>  (25) 

o    Yo 
Y 

<J>    =  -s  (26) 

o      o 
z 

u    +  v    =0  (27) 

o     o 
x      y 

These  show  the  flow  to  be  geostrophic,  hydrostatic,  and  hori- 
zontally non-divergent  to  a  zeroth  order  approximation.   The 
equations  to  first  order  in  Rossby  number  are: 

u        +uu        +vu        -   v,    -    3yv      =    -d),  (28) 

o,  o   o  o   o  1  J    o  Y] 

t  x  y  x 

v        +  u  v        +vv        +   u,    +    3yu      =   -<J>,  (29) 

o ,  o   o  o   o  1  2    o  ^  ly 

t  x  y  J 

♦l      ■   "Vo   -    s.l  (30) 

z  K 

ui  +  vi  +  r  <wips'z  =  °  (31) 

x  y  s  z 

[6    \ 

s         +us         +vs         +      -H    w,    =    0  (32) 

o.  o   o  o   o  \  u  /       1 

t  x  y 

To  eliminate  the  first  order  pressure  terms  from  (28) 
and  (29)  the  equations  are  cross-differentiated  and  sub- 
tracted, resulting  in: 

(v      -u      )       +    u    (v      -u      )    +   v     (v      -u       )    +    (v      -u      )  (u      +v      ) 
o      "  o      t  oo         o  oo         o  o        oo        o 

xyu  xy  xy.xyxy 

+     (u,    +vn     )    +    py(u      +v      )    +    3v      -    0  (33) 

1        1  J      o        o  o 

x        y  x        y 

Then,  substitution  from  (27)  and  (31)  into  (33)  results  in: 
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(v   -u   )   +  u  (v   -u   )  +  v  (v   -u   ) 
o   o  t    oo   o      o  o   o 
x    y  u         x    y  x    y 


-=-  (w,  p    +  3v  =  0 
p     V  s'  v  o 

s       z 


(34) 


Equation  (32)  is  rewritten  in  the  following  form: 


yp 


(wlps}    +      6 


s 


(s        +us        +vs      )=0 
/  o,  o   o  o   o 

p  t  x  y 


PP. 


(35) 


Substitution  of  (26)  into  (35)  (noting  that  -? 7^  f  (x,y,t)) 

P 
gives : 


(wlPs) 


yp    \  fyp     \  /yp 

—  d>     +  u  ■? —  6    +  v  -? — - 

o     o   ,     o  I  o     o      o  1 6     o  , 


(36 


p     z'  t 


P     z  X 


P     z'  y 


and  differentiating  this  with  respect  to  z  gives: 


(wlpsl  = 


/yps        \  fyps        1 

7 d)        +  U   7 d)        +  V 

0     o  /   .    o  o     o  /      o 
P     z'  zt     l  P     z'  zx 


yp. 


6   Yo 
p    z  zy 


+  Uo    6—^0     +  Vo 
z  '  p     z'  x      : 


/ypQ 

o    o 
p    .z' 


(37) 


y 


which  can  then  be  substituted  directly  into  (34).   In  writing 
this  it  is  helpful  to  define  all  of  the  velocities  in  terms 
of  the  pressure,  assuming  a  geostrophic  relationship  (this 
assumption  filters  out  gravity  waves) . 


u   =  -d> 
o     ro 


y 


(v   -u   )  =  d>    +  d) 
°x   °y      XX     ^ 


o    Yo 


X 


(38) 


From  these  it  can  be  seen  that  the  last  two  terms  in  (37) 
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identically  cancel.   Then,  substitution  of  (37)  and  (38) 
into  (34)  and  dropping  all  zero  subscripts  gives: 


^t     x^y 


<J>  q   =  0 
Yynx 


(39) 


where : 


'yp, 


q  =  A    +  (J)    +  By  -  — 
M    vxx    Yyy    M^    ps 


(40) 


Equation  (39)  is  the  quasi-geostrophic  vorticity  equa- 
tion and  (40)  is  the  potential  vorticity.   The  first  two 
terms  on  the  right  hand  side  of  (40)  constitute  the  relative 
vorticity.   The  third  term  is  the  contribution  due  to  the 
variation  in  the  Coriolis  parameter,  and  the  last  term  repre- 
sents the  vorticity  contributed  by  vertical  stretching  of  the 
water  column.   Equation  (39)  shows  that  time-dependent  motion 
will  result  when  the  motion  is  such  as  to  create  a  change  in 
potential  vorticity.   For  the  initial  problem  studied  here, 
the  Coriolis  parameter,   6,  was  assumed  to  be  non-zero  and 
the  resultant  motion  investigated.   It  should  be  noted  that 
as  the  North  pole  is  approached,  3  approaches  zero,  but  that 
motion  could  still  be  driven  by  bathymetric  variations. 

Assuming  a  separable  wave  solution  of  the  form: 

$    =  Z(z)ei(kX+ly"at)  (41) 

leads  to  the  following  Sturm-Liouville  problem  for  the  verti- 
cal dependence: 


'^  Z  I   +  [k2  +  I2  +  & 

6P     z/z     L 


Z  =  0 


(42) 
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Substituting  from  eq .  (7)  for  p  and  6  ,  using  the  fact  that 


the  Brunt-Vaisala  frequency,  N,  is  given  by 
converting  back  to  a  dimensional  form  results  in: 
,    2 


g   Hs 


pq  8z 


and 


o 


s   u 


N' 


k2  +   2  +  Bk 

a 


Z  =  0 


(43) 


B.   EQUATION  USING  ARCTIC  DFNSITY  PROFILE 
1 .   Basic  State  Description 

In  this  initial  formulation  zero  mean  flow  was 
assumed,  thus  the  initial  density  distribution  was  assumed 
to  be  a  function  of  depth  only.   The  Arctic  under-ice  density 
profile  was  approximated  by  a  homogeneous  upper  layer  from 
the  surface  to  z  =  -h,  below  which  density  increased  follow- 
ing an  exponential  curve.   In  this  lower  layer  density 
increases  such  that  it  assymptotically  approaches  a  constant, 
p , ,  at  deeper  depths.   This  distribution  is  similar  to  that 
proposed  by  Reid  [1948] ,  in  his  equatorial  model,  and  is 
given  analytically  by: 


Ks    "n 


ps  =  pd  "  (pd"ph)e 


(0  >  z  >  -h) 


(44) 


l/c(z+h) 


(z  <  -h) 


which  is  illustrated  in  Figure  8  for  various  values  of  c  and 
h  (referred  to  hereafter  as  layer  depth  (h)  and  e-fold  depth 
(c) ) .   Actual  density  values  taken  during  the  collection  of 
velocity  data  are  also  shown  on  Figure  8  for  comparison.   As 
can  be  seen,  the  best  data  fit  is  with  an  e-fold  depth  of 
between  100  and  120  meters  and  a  layer  depth  of  30  meters. 
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Figure  8.   Depth  Distribution  of  Analytic  Density 
Function  Shov.'ing  Relation  to  Data. 
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2 .   Use  of  Actual  Density  Distribution  in  Equations 

The  density  distribution  given  by  eq.  (44)  was  sub- 
stituted into  eq.  (43)  to  produce  a  vertical  normal  mode 
equation  which  would  apply  to  this  specific  density  structure. 
To  show  this,  eq .  (43)  is  first  expanded  in  terms  of  p   and 
its  derivatives : 


o 


-^—   Z         +    2Z       - 
p  zz  z 


L.     S 


p    p 
s*s 


zz 


Ups   ) 

z 


ZJ 


k2.+l2  +  ^ 

a 


Z  =  0     (45) 


(46) 


The  derivatives  of  p   are  found  from  eq .  (4  4)  as 

"^d"^   l/c(z+h) 
p_   =  e 


(pd  ph}   l/c(z+h) 

=  2 e 

zz       c 


which  can  then  be  substituted  into  eq.  (44)  to  yield,  after 
some  rearrangement  of  terms : 


Z 


zz 


"2(pd-ph}   l/c(z+h)  M   1 

e         +  — 

p  C  CJ 


f   c 
o 


2~  (pd-ph' 


k2  +  I2  +  to 


el/c(z+h)z 


(47) 


This  can  be  simplified  to: 


zz 


Ael/c(z+h)  + 


z   +  XeVc(z+h)    .  Q 


(48) 


where 


A  = 


2(pfl-Ph) 


P  c 
s 


(49) 
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B  = 


(49) 


2 

f   cp 
o   Ks 


(pd-ph} 


2     2 
k  +  1  + 


Pk 


The  solutions  to  eq.  (48)  will  then  be  a  set  of  eigen- 
functions  describing  the  vertical  distribution  of  the  pressure 
function,  $ .   It  is  noted  that  this  equation  only  applies  to 
the  region  below  the  layer  depth.   In  the  upper  homogeneous 
region  horizontal  velocities  will  be  constant  since  there  is 
no  stratification.   Vertical  boundary  conditions  on  the 
system  are  that  w,  the  vertical  velocity,  is  equal  to  zero  at 
both  the  surface  (ice)  and  bottom  boundaries.   In  addition, 
velocities  are  assumed  to  be  continuous  across  the  internal 
boundary  at  z  =  -h.   If  these  boundary  conditions  are  met 
and  a  solution  to  eq.  (48)  exists  there  will  be  solutions 
for  different  discrete  values  of  A,  the  eigenvalues  for  this 
equation. 

A  value  of  either  k  or  a  may  then  be  assumed  and  used 
with  the  eigenvalues  to  solve  for  the  corresponding  values 
of  a  or  k  for  each  mode.   In  order  to  actually  solve  for  both 
k  and  a  some  initial  distribution  for  <j)  must  be  known  and 
matched.   This  is  the  initial  value  problem,  based  on  the 
perturbation  at  a  lead,  which  is  the  next  phase  of  the  study 
and  is  discussed  in  Section  IV. 

The  problem  being  discussed  here  was  formulated  with 

the  assumption  of  a  long  narrow  lead  aligned  in  a  North-South 

direction.   Thus  x,  «  y. ,  where  x,  and  y,  are  the  Fast-West 

1    J 1  1      Jl 
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and  North-South  dimensions  of  the  lead,  respectively.   Based 
on  this  it  was  then  assumed  that  no  changes  occur  in  the  y 
direction  (9/8y  =  0)  which  in  turn  means  that  1  is  assumed 
to  be  zero  in  eqs .  (41)  and  (43)  and  that  u  =-8cJ)/3y  =  0. 
The  equations  thus  become  two-dimensional.   A  horizontal 
boundary  condition  at  time  zero  of  v  =  0  at  the  origin, 
which  is  placed  in  the  center  of  the  lead,  is  then  imposed. 
The  horizontal  velocity  can  be  written  as: 

d<b         .  i  _  ,  .  i  (kx-ot) 
v    =    rr1-  =    ikZ  (z)e 
8x 

or : 

v  =  kZ  (z)  [C, cos (kx-at)  +  C  sin (kx-at)] 

For  v(0,0)  =  0,  C,  must  equal  zero  or: 

v  =  kZ (z)sin  (kx-at)  (50) 

This  means  that  c|>  can  be  written  as: 

(J)  =  Z  (z)  cos  (kx-at) 
From  eq.  (36)  w  can  be  written  as: 

yps 

w  =  r Z  sin (kx-at) 

0     z 

p 

which  becomes : 

f 
w  =  — y  Z  sin (kx-at)  (52) 

IT   z 

when  the  terms  are  dimensionalized  and  proper  substitutions 

made  for  y  and  6  .   Assuming  a  value  for  k  is  known,  an 
expression  for  a  can  then  be  derived  from  eq .  (49) : 


35 


0  =   —2 (53) 

f   cp  X      „ 

g(Pd-Ph) 

This  is  similar  to  the  equation  for  the  frequency  of  a  baro- 
tropic  Rossby  wave,  a,  =  $/k ,    v;ith  the  additional  term  in 
the  denominator  being  due  to  the  stratification. 
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III.   RESULTS  OF  WAVE  SOLUTION  AFP  COMPARISON 
WITH  ACTUAL  DATA 

A.   NUMERICAL  SOLUTION  OF  THE  VERTICAL  MODE  EQUATION 

The  vertical  mode  equation  (48)  was  rewritten  using  a 
normalized  depth  variable,  0  >_  z '  >_  -1.0,  as: 

I    -  D(AeDB(z'+h/D)+B)Z-  -X'eDB(z'+h/D»Z=0     (54) 
zz  z 

2 

where  X1  =  D  X  and  D  =  depth  of  water. 

Equation  (54)  was  solved  numerically  on  the  Naval  Post- 
graduate School  IBM/360  computer.   The  solution  was  found 
using  a  Milne  predictor-corrector  method,  according  to  an 
algorithm  described  by  Nielsen  [1964].   The  depth  range  was 
divided  into  2000  increments  and  then  the  Milne  method  used 
from  the  bottom,  z1  =  -1.0,  upward  to  the  normalized  layer 
depth,  z1  =  -  h/D.   The  Milne  method  section  of  the  numerical 
solution  was  then  tested  by  using  it  to  solve  a  zero  order 
Bessel  equation  which  had  a  known  analytic  solution.   Results 
of  this  test  were  satisfactorily  checked  with  a  table  of 
Bessel  functions. 

In  the  relatively  thin  homogeneous  region  above  the  layer 

depth,  the  slope  of  the  vertical  velocity  was. held  at  its 

layer  depth  value,  as  was  the  value  of  the  vertical  structure 

function,  Z.   The  output  of  the  numerical  scheme  was  a  set  of 

2000  values  of  Z  and  Z   from  z1  =  -1.0  to  z'  =  0,  for  a  given 

z 

X*.   The  vertical  distribution  of  horizontal  velocity  is 
proportional  to  Z  as  was  shown  by  eq .  (50).   From  eq .  (52)  it 
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can  be  seen  that  the  vertical  dependence  of  vertical  velocity, 

2 
w,  is  proportional  to  7,    /N  .   The  numerical  scheme  could  thus 

be  used  to  give  the  vertical  dependence  of  both  horizontal 

and  vertical  velocities. 

The  bottom  boundary  condition,  w  =  0  at  z'  =  -1.0,  was 

met  by  setting  Z  (-1.0)  =  0.   This  was  valid  since  w  is 

2  2 

proportional  to  Z  /N   and  since  N  (-1.0),  although  quite 

small,  is  not  identically  zero.   With  the  bottom  boundary 

condition  thus  fixed,  the  condition  w  =  0  at  z '  =  0  could 

only  be  met  for  certain  values  of  X  '  .   These  values  are  the 

eigenvalues  corresponding  to  the  normal  modes  of  oscillation 

of  the  ocean  with  the  given  vertical  density  distribution. 

In  order  to  find  the  correct  eigenvalues  an  initial  X '  was 

assumed.   An  iterative  process  was  then  begun  in  which  the 

numerical  method  described  earlier  was  used  to  arrive  at  a 

value  of  w(0)  for  each  new  value  of  X ' .   The  initial  value 

chosen  for  X'  was  then  increased  in  steps  until  the  value 

of.  w(0)  changed  sign.   At  that  point  a  Regula-Falsi  method 

(eq.  (55)  below)  was  used  to  find  the  next  X'. 

X'   w  (0)  -  X'w   ,  (0) 
\  >    -  n~1    n n  n-1  ,  r  c-  x 

n+1       w  (0)  -  w   ,  (0)  KDD> 

n        n-1 

This  was  continued  until  the  value  of  w(0)  was  sufficiently 
close  to  zero.   The  value  of  X '  at  that  point  was  then  con- 
sidered an  eigenvalue.   The  process  was  continued  by  increas- 
ing the  value  of  X '  in  steps  until  the  next  eigenvalue  was 
found  in  the  same  manner. 
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B.   RESULTS  OF  THEORETICAL  DFVELOPMFNT 

Figures  9  and  10  show  the  vertical  distribution  of 
horizontal  and  vertical  velocities  which  result  from  the 
numerical  solution  of  eq.  (54)  for  a  given  density  distribu- 
tion (h=30,  c=120) .   The  first  four  normal  modes  are  shown 
in  each  figure.   Velocity  values  plotted  are  not  absolute 
values  but  are  relative  to  the  horizontal  velocity  at  the 
bottom  for  each  mode.   The  ratio  of  w(z)  to  V(D)  was  computed 
using  a  and  k  values  corresponding  to  a  wavelength  of  30  km. 
As  can  be  seen  on  these  graphs,  the  vertical  velocity  would 
be  insignificant  for  a  horizontal  velocity  on  the  order  of 
cm/sec  or  larger.   This  means  the  only  measurable  motion 
would  be  in  the  horizontal  direction. 

Figure  11  shows  the  effect  of  a  change  in  the  density 
distribution  layer  depth  on  the  horizontal  velocity  distribu- 
tion.  The  second  mode  only  is  shown  since  effects  on  other 
modes  are  similar.   The  values  chosen,  20  and  50  meters, 
represent  the  normal  seasonal  variation  in  the  thickness  of 
the  under-ice  layer.   As  can  be  seen,  very  little  change  in 
the  velocity  distribution  results. 

Figure  12  shows  the  effect  on  the  second  mode  of  a  change 
in  the  curvature  (e-fold  depth,  c)  of  the  below  layer  density 
distribution.   This  shows  that  as  the  e-fold  depth  increases 
(decrease  in  curvature) ,  the  velocity  peak  is  less  marked 
and  occurs  at  a  deeper  depth. 

Frequencies  were  solved  for  by  using  eq .  (53)  with  an 
assumed  horizontal  wavelength  of  30  km.   Results  were  on  the 
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Figure  9.   Vertical  Structure  of  Theoretical  Horizontal 
Velocity  Patterns  for  Modes;  1  through  4. 
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Figure  10. 


Vertical  Structure  of  Theoretical  Vertical 
Velocity  Patterns  for  Modes  1  through  4. 
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Figure  11.   Effect  of  Layer  Depth  Change  on  Vertical  Structure 
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Figure  12. 


Effect  of  Varying  e-fold  Depth  in  Density 
Distribution  on  Vertical  Structure. 


—  8     —  1 
order  of  10    sec    or  less,  which  corresponds  to  periods  of 

tens  of  years.   These  values  indicate  the  current  patterns 

are  essentially  stationary.   As  can  he  seen  from  eq .  (53) 

either  a  decrease  in  wavelength  (increase  in  k)  or  an  increase 

in  mode  (higher  A)  will  cause  the  frequency  to  decrease. 

Thus  the  limiting,  highest  possible  frequency  for  a  given 

density  distribution  and  horizontal  wavenumber  is  for  A  =  0, 

a  barotropic  Rossby  wave  mode,  and  is  given  by: 

b    k 

For  a  wavelength  of  30  km  this  corresponds  to  a  period  of 
about  10  years. 

The  horizontal  distribution  of  velocities  was  arrived  at 
by  assuming  a  length  of  15  km  to  the  first  zero  of  the 
sin(kx-ot)  function  of  eq.  (50)  and  forming  the  product  of 
this  trigonometric  function  and  the  Z (z)  function  at  a  con- 
stant depth.   This  is  shown  in  Figure  13. 

C.   COMPARISON  WITH  OBSERVED  DATA 

A  comparison  of  the  theoretical  normal  modes  plotted  in 
Figure  9  with  the  observed  data  in  Figure  2  shows  that  the 
second  normal  mode  most  nearly  matches  the  measured  vertical 
distribution  of  horizontal  velocity.   The  only  real  discre- 
pancy between  the  two  is  the  relative  maximum  in  the  theo- 
retical distribution  at  the  surface.   This  appears  to  be  an 
inherent  characteristic  of  the  theoretical  solution  and  is 
seen  in  each  of  the  normal  modes.   One  possible  explanation 
is  that  these  theoretical  maxima  were  caused  by  treatment  of 
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Figure  13.   Horizontal  Distribution  of  Theoretical 
Horizontal  Velocities. 
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the  under-ice  surface  as  flat  and  frictionless  in  the  develop- 
ment.  The  ice  actually  has  a  rough  surface  capable  of  dis- 
sipating some  energy  of  relatively  slow  moving  currents. 

A  second  explanation  is  that  more  than  one  mode  exists 
with  a  progressively  decreasing  amount  of  energy  in  succes- 
sively higher  modes.   Thus  if  v  (z)  represents  the  horizontal 
velocity  distribution  for  a  given  mode,  then: 

oc 

v  (z)  =  S   a  v  (z  )  , 
,   n  n 
n=l 

where  a   >  a  ,,,  would  describe  the  actual  velocity  distri- 

n     n+1  J 

bution.   Assuming  further,  as  a  simplistic  first  approach, 
that  the  only  significant  a   are  for  n=l  and  n=2  (i.e.,  only 
the  first  two  modes  contain  significant  energy)  the  linear 
combination  of  modes  1  and  2  was  plotted  in  Figure  14.   The 
a,  and  a„  were  sized  to  give  zero  velocity  at  the  surface, 
meaning  a~  was  less  than  a, .   Figure  14  can  be  seen  to  be 
closer  to  the  observed  data  than  the  second  mode  alone. 
Obviously  a  more  rigorous  solution  for  the  coefficients  could 
be  used  to  get  an  even  closer  fit  to  the  data.   This,  however 
would  not  be  physically  meaningful.   The  purpose  of  this 
initial  investigation  was  merely  to  show  that  the  theoretical 
development  would  allow  motion  similar  to  what  has  been 
observed.   This  has  been  shown.   The  next  phase  of  the  inves- 
tigation, which  will  be  discussed  in  Section  IV  of  this 
thesis,  is  to  study  the  initial  value  problem  of  how  this 
specific  ocean  system  responds  to  the  initial  perturbations 
in  its  density  structure  caused  by  a  newly  forming  lead 
(i.e.,  which  normal  modes,  if  any,  are  excited). 
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Figure  14.   Linear  Sum  of  Fodes  1  and  ?    for 
h  =  30,  c  -  120. 


In  order  to  compare  real  and  theoretical  distributions 
more  closely,  a  vertical  cross-section  through  the  origin 
was  drawn  showing  theoretical  horizontal  velocity  values, 
Figure  15.   This  was  done  for  the  second  normal  mode  for  a 
density  distribution  having  h=30  and  c=120,  and  for  the 
simple  horizontal  distribution  of  Figure  13.   Thus  Figure  15 
could  be  compared  with  Figure  4,  which  is  drawn  for  real 
current  data  having  a  similar  density  distribution.   A 
similarity  can  be  seen  in  the  general  flow  direction  and 
in  the  shape  of  the  isolines.   Apart  from  the  relative 
smoothness  of  the  theoretical  plot,  the  only  significant 
difference  is  in  the  upper  homogeneous  layer,  which  has  been 
discussed  previously.   Figure  16  is  the  same  cross-section 
drawn  for  the  summation  of  the  first  two  normal  modes.   As 
can  be  seen,  the  cell  structure  is  closer  to  the  observed 
data  of  Figure  6  than  was  the  second  mode  alone. 
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Figure  15.   Vertical  Cross-section  of  Horizontal 
Velocities,  Mode  2. 
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Figure  16.   Vertical  Cross-section  of  Horizontal 
Velocities,  Sum  of  Modes  1  and  2. 
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IV.   DEVELOPMENT  OF  THE  INITIAL  VALUE  PROBLEM 

A.   GENERAL  METHOD 

The  problem  addressed  here  is  whether  the  density  field 
perturbations  which  occur  at  a  newly  forming  Arctic  lead  can 
force  the  theoretically  developed  oscillations  to  occur  and, 
if  so,  whether  the  forced  oscillations  are  consistent  with 
observed  data. 

The  solution  to  this  problem  lies  in  formulating  expli- 
citly the  initial  effect  of  the  lead-induced  perturbations 
on  the  pressure  field.   This  effect  can  then  be  used  as  an 
initial  condition  on  the  pressure  function,  cf> ,  in  order  to 
fully  solve  the  Sturm-Liouville  problem  developed  in  Sections 
II  and  III.   It  was  shown  that  the  general  solution  to  that 
problem  was  a  set  of  eigenfunctions  which  gave  the  initial 
response  of  the  Arctic  Ocean  to  a  baroclinic  perturbation. 
The  linear  combination  of  these  functions  can  be  written  as: 

oc     <x 

d>(x,z,t)  =   2  '   S   A  B  Z  (z)cos(k  x-a   t)       (56) 
Y  .,       n  m  n         m   ran 

n=l  m=l 

so    that   at    time    zero: 

oo  o: 

<f>     (x,z)    =      2         £       A    B    Z     (z)cos(k    x)  (57) 

Yo  ,         ,       n  m   n  m 

n=l   m=l 

It  must  be  noted  that  the  summation  of  solutions  indicated  by 
eq .  (56)  is  not  necessarily  a  solution  itself.  It  would  only 
constitute  a  solution  if  the  initial  equations  were  actually 
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linear,  which  they  are  not.   It  appears  reasonable  to  assume 

however  that  the  non-linear  contribution  is  small  compared 

to  the  linear  so  that  this  summation  approximates  an  actual 

solution.   A  method  will  now  be  shown  whereby  the  values  of 

A  ,B  ,k   and  a    can  be  found  based  on  the  initial  specifica- 
n   m   m       mn  L 

tion  of  the  <b      field, 
o 

B.   SPECIFICATION  OF  INITIAL  CONDITIONS 

The  pei"turbations  at  a  lead  during  formation  and  refreez- 
ing  were  described  as  small  downward  velocities  and  an 
increase  in  density  under  the  lead  due  to  salt  rejection  and 
small  downward  velocities  due  to  shear  processes.   Consider- 
ing the  increase  in  density  below  the  lead  as  the  dominant 
process,  the  initial  density  perturbation  will  result  in  a 
deviation  from  the  basic  state  as  shown  in  Figure  17;  where 
z   shown  is  the  depth  of  penetration  of  the  density  pertur- 
bation.  Assuming  that  this  takes  place  only  directly  under 
the  lead,  then  a  vertical  cross-section  of  initial  density 
would  appear  as  shown  in  Figure  18.   Using  the  hydrostatic 
assumption,  this  will  result  in  an  initial  overpressure  con- 
dition (perturbation  pressure  p'  =  $    )  as  shown  in  Figure  19 

for  -x   <  x  <  x   and  zero  elsewhere.   This  can  then  be  used 
p  -   -  p 

as  an  initial  condition  to  which  eq.  (57)  can  be  fitted  in 
order  to  solve  for  the  appropriate  wavenumbers ,  frequencies, 


v   +-   |z  -h 

bat  '  p 


and  constants.   A  further  assumption  that  '  p  '  '  <<  1,  will 

h 

be  made,  thus  z   ~  -h  can  be  used  in  describing  the  initial 

P 

distribution  analytically.   (J)  (x,z)  can  thus  be  written: 
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Figure  17.   Perturbation  Density  Due  to  Salt 
Rejection,  at  Tine  Zero. 
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Figure  18.   Vertical  Cross-section  of  Perturbation 
Density  at  Time  Zero. 
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Figure  19.   Perturbation  Pressure 
at  Time  Zero. 
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X    *^  X 

♦0(x,z)  =  -  ~     I  '       P  (58) 

0  <  z  <  z 
-   -  p 

X  I   <  X 

=  1    I     -  p 

z   >  z 

p; 

=0      ( Ixl  >  x  ) 

1  '   p' 


C.   SOLUTION  OF  THE  PROBLEM 

The  initial  condition,  cf>  (x,z),  is  separable  and  can  be 
written  as  <j)  (x,z)  =  f  (x)g(z)  where: 

f(x)  =  0    (|x|  >  x  )  (59) 

=1    ( Ixl  <  x  ) 


g(z)  =  -z/h    (z  >  -h)  (60) 

=1      (z  <  -h) 

Thus,  coefficients  can  be  solved  for  separately  from: 


f  (x)  =  Z       B  cos  (k  x) 
,  m  m 

m=l 

(61) 

oo 

g  ( z )  =   E   A  Z  ( z ) 
t   n  n 
n=l 

Considering  B   first,  these  are  the  coefficients  for  a 
Fourier  cosine  series,  which  on  the  interval  (0,L)  are  given 

by: 

x 
-i    P  i   L  x 

B^  =  -  /    f  (x)dx  +   =-  f      f  (x)dx  =  -£■  (62) 

°    L  0  L  x  L 

x  p 

Bm  -  |  /  P  f  (x)cos(~^)dx   (111=1,2,3...) 

=  ~   sin(kmxp)     (m=l,2,3...)  (63) 
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where  k   =  miT/L  is  a  distance  sufficiently  far  removed  from 

m  J 

the  lead  so  that  4>  can  be  assumed  to  have  a  zero  value. 

The  values  of  B   correspond  to  the  zeroth  mode,  which  is 

o       ^ 

a  barotropic  mode  that  may  exist  but  is  not  of  interest  in 
the  present  discussion.   Figure  20  is  a  plot  of  the  approx- 
imation to  f (x)  formed  by  evaluating  eg.  (61)  for  m  =  1 
through  m  =  40.   It  can  be  seen  to  approximate  f (x)  fairly 
closely.   Figure  21  is  a  plot  of  the  horizontal  distribution 
of  v  for  these  same  modes. 

For  each  k   thus  found,  and  each  eigenvalue  corresponding 

to  the  Z  (z)  vertical  structure  functions,  there  can  be  found 
n 

a  value  of  a 


mn 

The  values  of  A   are  given  by: 

n      ^       J 


A   = 
n 


-D 
/   g(z)Z  (z)dz 
0 


-D 


(n=l , 2 , 3 . . . ) 


(64) 


2 


/.    [Z  (z)]"dz 
0      n 


for  a  general  series  expansion,  assuming  the  Z  (z)  are 

orthogonal  functions  with  respect  to  a  weighting  function 

of  one.   The  denominator  can  be  normalized  for  all  modes  by 

finding  a  factor  C   such  that: 

n 

"D  2 

/  [C    Z     (z)rdz    =    1 

o  n  n 


or : 


n 


"D  2      ' 

/     (?,  u)rdz 

-  o 


(65) 


Thus   A      becomes: 
n 
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Figure  20.   Fourier  Series  Fxpansion  of  Initial 
Horizontal  Pressure  Distribution. 
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Figure  21.   Fourier  Series  Fxpansion  of  Horizontal 
Velocity  Distribution  in  ,x-direction . 
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-D 

A   =  /    g(z)C  Z  (z)dz    (n=l,2,3...)  (66) 

n       -,     n  n 

for  the  normalized  distributions.   The  computation  of  A   can 
be  further  simplified  by  writing  g(z)  as  the  sum  of  an 
average  and  a  variable  part.   Since  D  »  h,  the  average  part 
is  essentially: 

g    (z)  -  1. 

^avg 

The  only  variable  part  is  then: 

g.ra„(z)  =  1  +  z/h    (z  >  -h)  (67) 

Thus,  eq.  (64)  can  be  written: 

-D  -h 

A      =   f        C  Z  (z)dz  +  /    (l+z/h)C  Z  (z)dz       (68) 
n    '     nn        '        '  '  n  n 

But,  the  first  integral  in  eq .  (68)  is  equal  to  zero  based  on 
the  original  continuity  relationship  used  in  deriving  Z(z). 

From  Figure  9  it  can  be  seen  that  each  Z  (z)  is  constant  for 

a  n 

depths  above  -h.   Thus  eq .  (68)  finally  becomes: 

A   =  -  h/2  [a  C  ]    (n=l,2,3...)  (69) 

n  n  n 

These  values  were  solved  for  using  the  IBM/360  computer  and 

were  then  used  to  form  the  series  summation  given  by  eq.  (61) 

for  the  first  six  modes.   The  results  are  shown  in  Figure  22. 

The  approximation  is  not  as  close  as  the  horizontal  case  due 

to  using  only  six  modes.   This  was  a  limit  imposed  by  the 

numerical  method. 

A  complete  set  of  values  for  A  ,  B  ,  k  ,  and  a    is  then 
^  n   m   m       mn 

available.   Values  for  the  distribution  having  c=120,  h=30, 
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Figure  22.   Series  Expansion  of  Vertical 
Structure  Function. 
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are  given  in  Tables  I,  II,  and  III.   Values  of  L=20  km  and 

r  =1  km  were  assumed. 
P 

The  solution  for  4>  given  by  eq .  (56)  will  then  begin 
with  a  distribution  as  given  by  eq .  (57),  which  is  the  ini- 
tial pressure  distribution.   Due  to  the  extremely  small 
values  of  a  there  will  be  no  significant  movement  with  time 

of  the  overall  wave  system.   Phase  velocities  (a/k)  are  on 

_2 

the  order  of  10    cm/sec.   The  distribution  of  horizontal 

velocities  however,  is  proportional  to: 

oo      oo 

v(x,z,t)  =   E    7,      A  B  k  Z  (z)sin(k  x-c   t)     (70) 
'  '  '      ,    ,nmmn        m   mn 
n=l  m=l 

The  initial  horizontal  and  vertical  distributions  of  v(x,z,t) 
were  shown  in  Figure  21  and  Figure  22. 

TABLF  I 

VALUFS  OF  o   (sec~1xl0~9)  AND  T   (yrs)  FOR 
mn  mn  J. 

SELECTED  VERTICAL  AND  HORIZONTAL  MODES 


m 


10 


1 

9.8  3 

5.84 

3.11 

0 

mn 

20.0 

33.9 

63.4 

T 

mn 

n      3 

2.58 

4.38 

2.85 

a 
mn 

76.9 

45.2 

69.3 

T 
mn 

5 

1.11 

3.0  2 

2.49 

a 
mn 

178. 

65.6 

77.0 

T 

mn 

60 


TABLE  II 


VALUES  OF  k   AND  B   FOR  HORIZONTAL 
m      m 

MODES  1  THROUGH  10 


m 

_  1    _  f- 
km(cm   xlO   ) 

Bm 

1 

1.56 

.0995 

2 

3.14 

.0983 

3 

4.70 

•  .0963 

4 

6.28 

.0935 

5    , 

7.85 

.0900 

6 

9.40 

.0858 

7 

10.1 

.0810 

8 

12.5 

.0756 

9 

14.1 

.0699 

10 

15.6 

.0637 

TABLE  III 


VALUES  OF  X   AND  A   FOR  VERTICAL 
n       n 

MODES  1  THROUGH  6 


n 


n 


n 


1 

+.0237 

9.18xl0-5 

2 

-.0218 

4.66xl0~4 

3 

+.0193 

1.14xl0~3 

4 

-.0169 

2.14xl0~3 

5 

+.0150 

3.47xl0~3 

6 

-.0134 

5.15xl0~3 
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V.   CONCLUSIONS 

The  stated  objective  of  this  study  was  to  model  certain 
observed  Arctic  Ocean  internal  current  patterns  through  an 
analytic  study  of  the  baroclinic  response  of  the  ocean  to 
lead  formation  and  refreezing.   The  problem  was  then  divided 
into  two  major  phases:   (1)  finding  the  baroclinic  internal 
normal  modes  of  the  Arctic  Ocean;  and  (2)  investigation  of 
the  forcing  of  these  modes  by  the  effects  of  lead  formation. 

The  normal  modes  of  response  found  as  a   result  of  the 
theory  developed  in  Section  II  showed  a  velocity  pattern 
which  was  essentially  stationary  with  no  significant  vertical 
velocity.   However,  horizontal  velocities  of  the  magnitude 
measured  in  the  Arctic  could  exist  within  these  stationary 
patterns.   It  was  noted  that  the  steepness  of  the  pycnocline 
(measured  by  the  parameter,  c)  had  a  direct  and  significant 
effect  on  the  depth  of  maximum  velocity,  whereas  the  effect 
of  a  variation  in  layer  depth  was  relatively  minor.   The 
apparent  consistency  of  the  cell  structure  seen  in  Figure  16, 
the  theoretical  vertical  cross-section  of  horizontal  veloci- 
ties ,  with  the  observed  data  lends  strength  to  the  plaus- 
ibility of  the  developed  theory. 

The  fact  that  a  simple  linear  combination  of  the  first 
two  theoretical  normal  modes,  Figure  14,  qualitatively 
matches  measured  velocity  profiles  is  of  particular  note. 
It  shows  the  physics  used  to  develop  the  theory  to  be 
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consistent  with  actual  data,  in  that  observed  current 
patterns  appear  to  be  part  of  the  set  of  allowed  solutions. 
It  has  not  yet  been  shown  conclusively  that  the  perturbations 
at  a  lead  will  force  these  modes. 

The  study  done  here  is  a  first  step  which  looked  at  the 
dynamics  of  a  single  Arctic  lead.   The  .theory  developed  has 
been  shown  to  be  dynamically  possible  but  the  study  has  not 
yet  looked  at  the  energy  involved  in  the  system.   It  is 
possible  that  an  interaction  between  leads  or  between  lead(s) 
and  the  mean  flow  is  needed  to  provide  the  energy  required 
for  the  observed  current  magnitudes. 

The  very  small  values  of  frequency,  corresponding  to 
periods  of  greater  than  20  years,  are  attributable  to  the 
magnitude  of  3  in  the  Arctic.   The  vorticity  equation  derived, 
eq.  (38),  is  a  non-linear  time  dependent  parabolic  equation. 
However,  the  results  found  in  the  initial  value  problem 
studied  in  Section  IV  of  this  thesis  indicate  that  the  time 
dependence  of  the  equation  is  not  significant  since  the  solu- 
tions were  essentially  stationary.   Under  the  same  geometry 
the  problem  could  be  reformulated  as  an  elliptic ,  linear 
equation  forced  with  time  dependent  boundary  conditions. 
An  advantage  of  this  type  formulation  would  be  the  ability 
to  write  a  separable  differential  equation  in  cylindrical 
coordinates ,  which  was  not  possible  using  the  development  of 
Section  II.   Motion  approximating  the  vortex-like  patterns 
observed  by  Coachman  and  Newton  [197  3]  could  then  be  more 
easily  described. 
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The  plane  wave  solution  for  velocities  developed  in  this 
thesis  is  a  result  of  the  formulation  of  the  problem  in 
rectangular  coordinates  with  an  infinite  length  lead.   This 
worked  well  in  developing  the  theory  and  resulted  in  a  real- 
istic vertical  section  of  velocity  structure.   However, 
actual  leads  are  of  finite  length  which  could  be  expected  to 
induce  curvature  in  the  velocity  patterns. 
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VI .   RECOMMENDATIONS  FOR  FURTHFR  RFSFARCII 

The  consistency  found  with  observed  data  warrants  the 
continuation  of  this  study  using  the  theory  developed  in  this 
thesis  as  a  basis.   The  following  recommendations  for 
further  research  are  made: 

1.  the  formulation  of  the  initial  value  "study  as  an 
elliptic  problem,  as  discussed  in  section  V; 

2.  the  use  of  vertical  stretching  caused  by  bathymetry, 
or  forced  convection  at  a  lead,  rather  than  the  3  effect  as 
the  driving  mechanism  in  the  developed  vorticity  equation; 

3.  the  investigation  of  the  question  of  whether  these 
velocity  patterns  have  any  influence  on  the  stability  of  the 
larger  scale  mean  flow,  including  the  effects  of  baroclinic 
instabilities ; 

4.  a  study  of  the  energetics  of  the  system. 


65 


BIBLIOGRAPHY 


Coachman,  L.  K. ,  "Production  of  Supercooled  Water  Durina  Sea 
Ice  Formation,"  Proc .  Symp .  Arctic  Heat  Budget  and  Atmos- 
pheric Circulation ,  The  Rand  Corporation,  RM-5233-NSF, 
497-529,  1966. 

Coachman,  L.  K.  and  J.  L.  Newton,  Personal  Communication, 
1973. 

Foster,  T .  D. ,  "Haline  Convection  Induced  by  the  Freezing  of 
Sea  Water,"  Jour.  Geophys .  Res.,  v.  73, 193 3- 1938,  March 
1968. 

Foster,  T.  D.,  "Haline  Convection  in  Polynyas  and  Leads," 
Jour.  Phys.  Ocean. ,  v.  2,  462-469,  October  1972. 

Gait,  J.  A.,  Current  Measurements  in  the  Canadian  Basin  of 

the  Arctic  Ocean,  University  of  Washington  Department  of 
Oceanography  Technical  Report  184,  1967. 

Gait,  J.  A. ,  "A  Numerical  Investigation  of  Arctic  Ocean 
Dynamics,"  Jour.  Phys.  Ocean.,  In  print,  October 
1973a. 

Gait,  J.  A.,  Unpublished  Manuscript,  1973b. 

Hildebrand,  F.  B. ,  Advanced  Calculus  for  Applications, 
Prentice-Hall,  Inc.,  1962. 

Kuo,  H.  L.,  "Dynamics  of  Quasigeostrophic  Flows  and  Insta- 
bility Theory,"  Advances  in  Applied  Mechanics,  v.  13, 
248-330,  Academic  Press,  1973. 

Nielsen,  Kaj  L.,  Methods  in  Numerical  Analysis,  MacMillan, 
1964. 

Pedlosky,  J.,  "The  Stability  of  Currents  in  the  Atmosphere 
and  Ocean:   Part  I,"  Jour.  Atmos .  Sci. ,  v.  21,  201-219, 
March  1964. 

Pedlosky,  J.,  "Geophysical  Fluid  Dynamics,"  Mathematical 
Problems  in  the  Geophysical  Sciences,  49-60,  American 
Mathematical  Society,  1971. 

Piascek,  S.  A.,  Numerical  Studies  of  Convection  Currents  at 
the  Ocean's  Surface  Caused  by  Fvaporation ,  Radiation  and 
Atmospheric  Cooling ,  Final  Technical  Report  for  ONR 
Contract  N0"0014-67-A-0242-003  ,  19  70. 


66 


Reid,  R.  O.,  "A  Model  of  the  Vertical  Structure  of  Mass  in 
Equatorial  Wind-driven  Currents  of  a  Baroclinic  Ocean," 
Jour.  Marine  Res.,  v.  7,  304-312,  November  1948. 

Schaus ,  R.  H. ,  A  Thermodynamic  Model  of  a  Central  Arctic  Open 
Lead,  M.S.  Thesis,  Naval  Postgraduate  School,  1971. 

Semtner,  A.  J.,  Personal  Communication,  1973. 


67 


INITIAL  DISTRIBUTION  LIST 

No.  Copies 


1.  Defense  Documentation  Center 
Cameron  Station 
Alexandria,  Virginia   22314 

2.  Oceanographer  of  the  Navy 
Hoffman  Building  No.  2 
2461  Eisenhower  Avenue 
Alexandria,  Virginia   22314 

3.  Office  of  Naval  Research 
Code  4  80  D 

Arlington,  Virginia  22217 

4 .  Naval  Postgraduate  School 
Code  58 

Monterey,  California   93940 

5.  Library,  Code  0  212 
Naval  Postgraduate  School 
Monterey,  California   93940 

6.  Dr.  R.  E.  Stevenson 
Scientific  Liaison  Office 
Scripps  Institution  of  Oceanography 
La  Jolla,  California   92037 

7.  Dr.  Jerry  A.  Gait 

NOAA  Pacific  Marine  Env.  Lab. 
University  of  Washington  WB-10 
Seattle,  Washington   98105 

8.  LCDR  Richard  I.  Itkin,  USN 
96  Bonaire  Circle 
Suffern,  New  York   10901 

9.  R.  K.  McGregor 

Director  Arctic  Program,  Code  415 
Office  of  Naval  Research 
Arlington,  Virginia   22217 

10.   AIDJEX 

Division  of  Marine  Resources 
University  of  Washington 
Seattle,  Washington   9  8105 


68 


11.  Assoc.  Professor  R.  T.  Williams 
Department  of  Meteorology 
Naval  Postgraduate  School 
Monterey,  California   93940 

12.  Professor  L.  K.  Coachman 
Department  of  Oceanography 
University  of  Washington,  WB-10 
Seattle,  Washington   98105 


6^ 


UNCLASSIFIED 


SECURITY   CLASSIFICATION   OF   THIS  PAGE   (Whan  Dale  Fntered) 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 

1.     REPORT  NUMBER 

2.  GOVT  ACCESSION  NO. 

3.     RECIPIENT'S  CATALOG   NUMBER 

4.     TITLE  (and  Subtitle) 

AN   ANALYTIC    MODEL    OF    OBSERVED    SMALL 
SCALE    BAROCLINIC    CURRENTS    IN    TEE 
ARCTIC    OCEAN 

5.     TYPE  OF  REPORT  &  PERIOD  COVERED 

Master's   Thesis; 
September    197  3 

6.     PERFORMING  ORG.    REPORT  NUMBER 

7.     AUTHORfs; 

Richard    Ivan    Itkin;    LCDR,    USN 

0.     CONTRACT  OR  GRANT  NUMBER(»J 

9.     PERFORMING  ORGANIZATION   NAME  AND  ADDRESS 

Naval   Postgraduate    School 
Monterey,    California      93940 

10.     PROGRAM   ELEMENT.  PROJECT,   TASK 
AREA  &   WORK  UNIT  NUMBERS 

II.     CONTROLLING  OFFICE  NAME   AND  ADDRESS 

Naval    Postgraduate    School 
Monterey,    California      93940 

12.     REPORT  DATE 

September    19  7  3 

13.     NUMBER  OF   PAGES 

71 

H.     MONITORING  AGENCY  NAME  A    ADDRESSf//  different  from  ConttolUn6  Office) 

IS.     SECURITY  CLASS,  (of  this  report) 

Unclassified 

15».     DECLASSIFI  CATION/ DOWN  GRADING 
SCHEDULE 

16.     DISTRIBUTION   ST  ATEMEN  T  (of  thte  Report) 

Approved    for   public    release;    distribution    unlimited. 

17.     DISTRIBUTION  STATEMENT  (of  tho  abstract  entered  In  Block  20,  If  different  from  Rnport) 

18.     SUPPLEMENTARY  NOTES 

19.     KEY  WORDS  (Continue  on  reverse  eldo  If  necessury  and  Identify  by  block  nu<r.ber) 

ARCTIC;    CURRFNTS;    BAROCLINIC;    QUASI-GEOSTROPHIC 

20.     ABSTRACT  (Continue  on  reverse  side  It  necessary  end  Identity  by  block  number) 

An   analytic   model    is    developed    of    baroclinic    currents 
measured   by    several    investigators    under    the   Arctic   Ocean 
ice    cover.      A    time-dependent,    Boussinesq,    quasi-geostrophic 
vorticity   equation    is    developed    to   describe    the    flow. 
Separable    normal   mode    solutions    are    found   which    show 
stationary    current   patterns    qualitatively    consistent  with 

DD  ,FJANRM73  1473 
(Page    1) 


EDITION  OF    1  NOV  65  IS  OBSOLETE 
S/N    0102-014-6601    I 

70 


UNCLASSIFIFD 


SECURITY   CLASSIFICATION   OF   KHIi  PAGE  (W„»n  Di.la  Entered) 


UNCLASSIFIED 


i'tCIJ'ilTY   CLASSIFICATION   OF   THIS   P  AGEOWion  Data  Enlorvd) 


observed  data.   An  initial  value  problem  is  then  formulated 
using  density  perturbations  at  a  nevzly  forming  lead  as  a 
forcing  function. 


DD     Form       1473  (BACK)  _.__ 

I  Jan  73  UNCLASSIFIED 

S/N       0102-014-6601  SECURITY   CLASSIFICATION   OF   THIS  PAGE(TW>»n  Dr.te  Entered) 

71 


'^H'15 


5 4593C 

,ode\   of 
e 
in 


^V* 


ltk»n 

&n  < 

obse    , .  r       cui  •  — 


*ed  -S*acurr< 


tbe 


Thesis 

145930 

175 

Itkin 

c.l 

An  analytic  model    of 

observed   small    scale 

baroclinic  currents    in 

the  Arctic  Ocean. 

thesl75 

An  analytic  model  of  nh 

mZmm^  fiscal 


3  2768  002  101 So 


